Quantum Hall plateau transition in the lowest Landau level of disordered graphene 
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We investigate, analytically and numerically, the effects of disorder on the density of states and 
on the localization properties of the relativistic two dimensional fermions in the lowest Landau 
level. Employing a supersymmetric technique, we calculate the exact density of states for the 
Cauchy (Lorentzian) distribution for various types of disorders. We use a numerical technique to 
establish the localization-delocalization (LD) transition in the lowest Landau level. For some types 
of disorder the LD transition is shown to belong to a different universality class, as compared to the 
corresponding nonrelativistic problem. The results are relevant to the integer quantum Hall plateau 
transitions observed in graphene. 



I. INTRODUCTION 

Recent experiments 1 - have unravelled a fascinating set 
of phenomena in atomically thin layer of hcxagonally ar- 
ranged carbon atoms known as graphene^ The quasi- 
particles of graphene are (2 + l)-dimensional massless 
Weyl ferrnionsi 3 -^ In the context of condensed matter 
physics their properties are strikingly different from non- 
relativistic fermions. And phenomena that are hard to 
realize for the relativistic case, such as the Klein paradox 
or the Zitterbewegung are accessible in graphene^ It is 
perhaps not an exaggeration to remark that many sub- 
tleties and a rich set of phenomenology are waiting to be 
discovered. 



A. Quantum Hall Effect 

A highlight has been the observation of an unconven- 
tional quantum Hall effectSiS and the corresponding the- 
oretical development ,2iI2iiiiI2ii2iI£ii5ii2iILi£ i n graphene 
the filling fractions are Vf = ±4(n + i) for magnetic field 
B < 9T, where n is an integer The factor of 4 comes 
from the two fold spin degeneracy and the two fold nodal 
degeneracy of the Landau levels. The Zeeman splitting is 
negligible compared to the cyclotron frequency and the 
disorder broadening of the Landau levels. The factor of 
half is due to a zero mode in the Landau level spectrum 
of Dirac fermions ^ 10 ' 11 

For stronger magnetic fields, 20T < B < 45T, plateaus 
appear at Vf = 0,±l,±2g, where q is an integer^ The 
plateaus at Vf — 0, ±1 can be explained by the lifting 
of both the spin and the nodal degeneracies in the low- 
est Landau level (LLL), but those at i>f = ±4, ±6, ... re- 
flect only the removal of spin degeneracy in higher Lan- 
dau levels. The removal of nodal degeneracy requires 
electron-electron interaction. Mechanisms suggested in- 
clude SU(4) ferromagnetis m 14 ' 15 , sublattice symmetry 
breaking due to short range interactions 1 ^ and the gen- 
eration of a mass gap by magnetic catalysis i 17 i 18 SU (4) 
quantum hall ferromagnetism predicts plateaus at all odd 



integer filling fractions. However, apart from v$ = ±1, 
the plateaus at Vf — ±3, ±5,.. have not yet been ob- 
served. 



B. Localization-delocalization transition 

The special quantization rules in graphene are ex- 
plained by the relativistic Landau levels, modified per- 
haps by interactions, but for the existence of Hall 
plateaus the Laughlin argument is necessary— Accord- 
ing to this argument the extended states at the center of 
a Landau band are separated by the localized states else- 
where. If the Fermi energy falls in the mobility gap, the 
plateaus are explained by a gauge invariance argument 
that is remarkably robust. The underlying phenomenon, 
therefore, is a localization-delocalization (LD) transition 
at the band cente n 19 ' 20 ' 21 ' 22 The conventional integer 
quantum Hall (IQH) plateau transition has been widely 
studied, and it is known that the localization length ex- 
ponent v « \ M Can we prove that the 

same argument applies to graphene, and, if so, does the 
LD transition belong to the same universality class? 



C. Disorder and Dirac fermions 

In the absence of a magnetic field, Dirac fermions in 
the presence of disorder have been widely studied in 
systems as varied as gapless semiconductors,— gapless 
superconductors ) 32 ' 33 and IQH plateau transitions ! 34 ' 35 
As compared to nonrelativistic fermions, the localiza- 
tion problem of Dirac fermions is richer because of a 
number of discrete symmetries. More specifically, if the 
disorder is particle-hole symmetric, for example a ran- 
dom gauge field, the LD transition takes place at zero 
energy and is reflected in the single particle density of 
states (DOS), in contrast to the conventional metal- 
insulator transition where the DOS is smooth through 
the LD transition. Surprisingly, there is a line of fixed 
points with continuously varying exponents depending 
on the disorder coupling constant 32 ' 33 ' 35 Some of the 
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unusual behavior of disordered Dirac fermions may be 
expected to realize in graphene. One such effect that 
has received considerable attention is the weak (anti)- 
localization phenomenon i 36 i 37 i 38 i 39 ' 40 However, relativis- 
tic Landau levels in the presence of disorder have not 
yet received much attention ! 11 ' 41 ' 42 ' 43 Here we provide a 
reasonably complete study of the possible effects. 

There is another important reason why LD transi- 
tions in the relativistic Landau level should be care- 
fully analyzed. In the conventional IQH effect, the spin- 
degenerate plateau transition corresponds to v ~ 4.6 
when it is assumed that the LD transition takes place 
at a single energy at the band cente n 44 ' 45 ' 46 ' 47 This has 
led to intense theoretical investigation of the LD tran- 
sition in the spin-degenerate Landau band i 48 ' 49 ' 50 ' 51 ' 52 
When spin-orbit scattering is included, the LD transi- 
tion is found to occur at two distinct energies, away from 
the band center. Scaling analysis about these distinct 
energies provide, once again, that v rs 7/3, as in the 
spin-polarized system. The scaling about a single energy 
at the band center leads to the effective exponent v ~ 4.6. 
One should anticipate a similar discrepancy between the 
spin and the nodal polarized IQH effect and the fourfold 
degenerate IQH effect in graphene. 



D. Graphene in the lowest Landau level 

For simplicity we shall concentrate on the spin polar- 
ized lowest Landau level (LLL) of graphene and analyze 
the LD transition both in the presence and in the ab- 
sence of nodal degeneracy. An interesting example of a 
controlled analytic calculation in the disordered Landau 
level problem is the DOS in the LLL. This was first com- 
puted exactly by Wegner— by examining the Euler trails 
of the impurity diagrams for the white noise disorder and 
was subsequently extended by Brezin et al£^ by using a 
supersymmctric (SUSY) technique. Here we also obtain 
some exact results for the DOS in the disordered rela- 
tivistic LLL using SUSY techniques. 

The most general model of disorder consists of a ran- 
dom potential, a random mass, a random gauge field 
and a random internode scattering; however, the random 
gauge field leaves the LLL unperturbed. After projection 
to the spin-polarized LLL, we study the following Hamil- 
tonian: 

3 

= mm + Vj (r)% , (1) 

3=0 

where Vo(r), V 3 {f) represent potential and mass disor- 
ders respectively and V\ (r) and V 2 (r) describe internode 
scattering effects. A mass m of the fermions have been 
included to study the effect of the removal of the nodal 
degeneracy. For simplicity we have omitted the constant 
Zeeman energy. The 2x2 matrix i]q is the identity matrix 
and rji, r\ 2 and 773 are the three Pauli matrices. 



E. Summary of results 

Because of a large number of cases involved, it is use- 
ful to summarize the results for the LD transition. Let 
9o, 93i 9i and 91 denote the widths of the Gaussian ran- 
dom distributions corresponding to the random poten- 
tial, random mass, and random internode scatterings, 
respectively. 

1. m = 



The list of possible cases are: 



f. 


.go 7^ and 53 = .91 = g 2 = 0. 




2. 


53^0, 


9o = 9i= 92 = 0. 




3. 


92 ¥= 0, 


.90 = .93 = .9i = 0. 




4. 


9i + 0, 


.9o = .93 = .92 = 0. 




5. 


.90 ± 0, 


.93 7^ and g x = g 2 = 0. 




6. 


50^0, 


g 2 7^ and g 3 = g x = 0. 




7. 


.90 ^ 0, 


.91 ^ and 53 = g 2 = 0. 




8. 


.93 ^ 0, 


.92 ^ 0, 50 = .9i = 0. 




9. 


.93 ± 0, 


.91 7^ 0, g = .92 = 0. 




10. 


.92 ± 0, 


.91 ^ and g = .93 = 0. 




11. 


.90 ? 0, 


.93 7^ 0, .92 7^ and g x = 





12. 


.90 ± 0, 


.93 7^ 0, .91 7^ and g 2 = 





13. 


50^0, 


92 ^ 0, gi ^ and g 3 = 





14. 


.93 ± 0, 


92 ^ 0, .91 ^ and .90 = 





15. 


.90 ? 0, 


.93 7^ 0, .92 7^ and .91 ^ 



In the cases (1) and (2), when disorder does not mix 
the two nodes, the LD transitions belong to the conven- 
tional IQH universality class with v sa 7/3. It is inter- 
esting to note that mass disorder produces LD transi- 
tion in the LLL, whereas for zero magnetic field random 
mass is known to be an irrelevant perturbation for the 
(2 + l)-dimensional Dirac fermions 1^ The Hamiltonians 
for (2), (3) and (4) involve only a single Pauli matrix at 
a time, related to each other by unitary transformations. 
Thus, (2), (3) and (4) are equivalent to each other and 
have v w 7/3. Because unitary transformations leave 
the identity matrix invariant, the same argument implies 
that (5), (6) and (7) are equivalent to each other and 
once again v « 7/3. 

The cases (8), (9) and (10) involve a pair of Pauli 
matrices and are equivalent to each other. In (8) the 
Hamiltonian has a discrete symmetry r\\Hr]\ = — H, of- 
ten called a particle-hole symmetry. The cases (9) and 
(10) have the same discrete symmetry with respect to 
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r]2 and 773. The case (10) has been analyzed by Hikami 
et al~ for a spin degenerate nonrelativistic LLL. When 
.9i = .92, DOS diverges at the band center and has two 
symmetrically located peaks away from it. The LD tran- 
sition takes place at these three distinct energies. Away 
from the band center the LD transition has the exponent 
v ~ 2.98 and the transition at the band center corre- 
sponds to a different exponent. If g\ ^ g 2 , the divergence 
of the DOS at the band center disappears, but the two 
symmetrically placed peaks away from the band center 
still exist. We find that the LD transition at these two 
energies hve continuously varying exponents depending 
on the ratio g 2 /gi- 

The cases (11), (12) and (13) are equivalent. The 
Hamiltonians in these cases are respectively the Hamil- 
tonians for the cases (8), (9) and (10), augmented by 
the identity matrix corresponding to the potential dis- 
order. Potential disorder breaks the discrete symmetry 
mentioned above, and there is no divergence of the DOS 
at the band center. The DOS is still peaked at two sym- 
metrically placed energies away from the band center. 
The LD transitions occur at energies away from the band 
center. If go is much smaller than the two remaining cou- 
pling constants, v follows trends similar to (8), (9), and 
(10). If go is comparable or larger, we find v ~ 7/3. 

In (14) all three Pauli matrices are present. The dis- 
crete symmetry of (8), (9) and (10) are absent, and the 
LD transitions take place at two symmetrically placed en- 
ergies away from the band center. When all the coupling 
constants are equal, the exponent v ~ 3.6. Depending 
on the relative strengths of the coupling constants the 
exponents vary continuously. If any particular coupling 
constant is significantly larger than the rest, v ~ 7/3. By 
adding go we obtain (15). If go is smaller than the rest, 
the situation is similar to (14). If go is larger than the 
rest, v ~ 7/3. 



2. m/0 

When m 7^ 0, the LD transitions occur at two symmet- 
rically placed energies about the band center, and these 
energies are greater than or equal to m. In the absence 
of internode scattering, the transitions occur at E = ±m 
and the exponent v ~ 7/3. If the strength of the intran- 
ode scattering is larger than m, the bands at ±m overlap 
and effectively correspond to the nodally degenerate case. 

If go = gz = and only one of the internode couplings 
is present, the DOS diverges at E — ±m with an ex- 
ponent of 0.5 and is identically zero for \E\ < m. The 
LD transitions occur at E = ±m and have a continuously 
varying exponent. When the disorder is strong compared 
to m, v ~ 7/3, and, in the opposite limit, v approaches 
unity. If we include small intranode scattering the situ- 
ation is similar. If the intranode scattering strength is 
greater than the internode scattering, v ~ 7/3. 

When go = gz — and both internode couplings are 
present, the DOS diverges at E = ±m with an exponent 



v r~j 0.47. However, the LD transitions occur at energies 
larger than \m\. We have analyzed a case where g\ = g 2 . 
The exponent varies continuously. If the internode scat- 
tering strength is larger than m, v ~ 3.8, and in the 
opposite limit v approaches unity. This behavior is sta- 
ble against intranode scattering if its strength is smaller 
than both m and the internode scattering. If intranode 
scattering strength is larger than the internode scatter- 
ing, v ~ 7/3. 

F. Roadmap 

Our paper is organized as follows: In Sec. II we de- 
scribe the Dirac fermion model. In Sec. Ill we describe 
various possible disorders and their forms when projected 
to the LLL. In sec IV we calculate the averaged density of 
states using supersymmetry. In the sections V, VI, and 
VII we describe the numerical studies of the LD transi- 
tion projected to the lowest Landau level. Section VIII is 
a brief concluding section. In the Appendix |A1 we provide 
some mathematical details of the density of states calcu- 
lation. In Appendix [B] we describe the recursive Green 
function technique used for numerical calculations and 
finally in Appendix [C] we outline the procedure of data 
collapse involved in the finite size scaling of the localiza- 
tion length. 

II. DIRAC FERMIONS AND LANDAU LEVELS 
OF GRAPHENE 

The low energy quasiparticles in graphene are well de- 
scribed by the Lorentz invariant form as the sum over two 
inequivalent nodes (the Fermi velocity vf ~ 10 6 m/s) 

H = -ifrvp J d 2 r§ CT (^D x + j 2 D y ) (2) 

where ^ a = Wj. 7 and the summation over spin a = ±1 
is understood. The four component Dirac spinor $ „ = 
{i'KA<y,i'4>KB<y,i^K l Ba,-'i^K l A<y), where the component 
ipKAa is constructed by superposing Bloch functions close 
to one of the two inequivalent nodes (K, K') of the Bril- 
louin zone, corresponding to one of the two sublatticcs 
(A, B) of the hexagonal graphene lattice. The notation 
D = (d — «f A) stands for the covariant derivative, A 
being the vector potential. The 7-matrices are defined 
by 7^ = (r 3 , iti, 1T2) <£> TI3] the Pauli matrix r operates on 
the two components corresponding to the sublattice in- 
dices, and the Pauli matrix rj operates on the components 
corresponding to the nodal indices. To be explicit: 

o As \ 1 (in \ 2 (ir 2 \ 

(3) 

To include a Zeeman term, we add 

H Z =E Z / 'dVj^Vg^'lV, (4) 
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where E z = g^sB is the Zeeman energy and Uj is a Pauli 
matrix operating on the spin indices. The Zeeman term 
breaks the SU (2) symmetry of the spin space down to 
U(l). The energy eigenvalues of the Hamiltonian opera- 
tor 

H = -ihvF^tfD* + i 2 D y ) + E z a 3 , (5) 
are well known: 



E nsa = s^2n\eB\Hv*/c-aE z ,n = 0,1,2,..., (6) 

where s = ±1 refer to the particle and the hole branches. 
In the presence of disorder Landau levels get broadened 
into a band, and the amount of broadening depends on 
the strength of the disorder. When the disorder is very 
strong, the half-width of the broadened band can be 
larger than E z , and experimentally this corresponds to 
the spin degeneracy of the Landau bands. In the spin 
degenerate situation, the observed filling factors is given 
b yI / f = 4(n+±)£I° 

The LLL wave function in the absence of the disorder 
in the symmetric gauge A = (—By/2,Bx/2,0) can be 
written as 



U(z, z) = e 




(7) 



where eB > 0. The functions fi(z) and f2(z) are holo- 
morphic functions of the complex coordinates z = x + iy; 
z — x — iy, and Ib = \J c/\eB\ is the magnetic length. 
Hence, in the zero mode, the first and the second node 
have nonzero amplitudes coming only from the sublat- 
tices A and B respectively. 

Two distinct onsite energies on the two sublattices cor- 
respond to a charge density modulation at the lattice 
scale. As a result, the particle and the hole branches 
acquire an energy gap. When linearized about the in- 
equivalent nodes, this energy gap appears as a parity 
preserving mass of the Dirac fermions. To be explicit, 
the linearized hamiltonian will have two new terms: the 
chemical potential term [(V4 + Vb)/2]#7°$ and the mass 
term [(Va — Vb)/2}^/^, where Va and Vb are the site en- 
ergies at the sublattices A and B. 

Although the non-interacting quasiparticles are mass- 
less in the absence of site modulation, they can acquire 
a parity conserving mass due to interaction effects. This 
spontaneous symmetry breaking is facilitated by the pres- 
ence of the magnetic field, a phenomenon known as "mag- 
netic catalysis" of chiral symmetry breakingi 56 i 57 The ef- 
fect has been argued to be the reason behind the quantum 
hall plateaus at Vf = 0,±1 observed in strong magnetic 
fields i 17 ' 18 Though it is beyond the scope of the present 
paper to consider electronic interactions, we will pay 
some attention to the noninteracting problem with a fi- 
nite mass. Our philosophy is to analyze the consequences 
of having a mass (possible in an interacting theory) on the 



LD transition. So, we shall include the term m i if'~f 0i $ in 
the effective Hamiltonian to examine the effect of mass. 
In the presence of such a mass term, the nodal degener- 
acy of Eq,s,o- = —cE z is removed and it splits into four 
levels £0,1, a = m — aE z and £0,— 1,0 = ~ &E Z . Each 
of these levels has the degeneracy if^. If the applied 
chemical potential is smaller than \E Z — m\, there will be 
a plateau at z/f = 0. If \E Z — m\ < < E z + m, = ±1 
plateaus will appear depending on the sign of ^i 16 i 56 Next 
possible values of quantized plateaus are V{ = ±2. The 
introduction of the mass term does not, however, lift the 
nodal degeneracy of the higher Landau levels, and the en- 
ergy levels i? n >i,a )<r = s^Jm 2 + 2n\eB\hvp / c — oE z has 

the degeneracy >^Sx, Therefore, when a mass is included, 
quantized plateaus appear at Vf — 0, ±1, ±2<7, where q is 
an integer. 



III. RANDOMNESS 

There are many sources of disorder in graphene : va- 
cancies, interstitials, substrate disorder and lattice dis- 
tortions due to dislocations. In principle there could 
also be random spin-orbit coupling. However, due to the 
small atomic mass of carbon, spin-orbit coupling is very 
weak compared to other energy scales. For simplicity, we 
shall primarily be interested in the spin polarized limit 
and ignore the random spin-orbit coupling. 

Point defects and substrate disorder can be described 
by introducing random site energies in the tight binding 
model. In the presence of substrate disorder there can 
also be a random modulation of the charge densities be- 
tween the two sublattices. These effects can be described 
by a random chemical potential Vb(r)^ , 7°4' and a ran- 
dom mass V3(r)^^> in the continuum limit. 

Because true long range crystalline order is not possi- 
ble in two dimensions at any finite temperature, topolog- 
ical defects, dislocations and disclinations will be present. 
Effects of these topological defects will result in random 
hopping amplitudes StAB and hence intranode as well as 
internode scattering. However, these scattering processes 
will take place between states on different sublattices. 
The follwing two bilinears, V 2 (r)^ , 7 3 'I' and Vi(r)*7 5 *, 
describe the internode scattering terms arising from ran- 
dom hopping. The two mutually anticommuting matri- 
ces, 
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I 

1 



,7 



-I 

1 



(8) 



also anticommute with 7 M ; I is the identity matrix. 

In the continuum limit, the most general impurity 
Hamiltonian is a 4 x 4 matrix: 

a -j a ^^ 2l(r) d 22 (t)J^ 

where Djj(r) are 2x2 matrices. Here D±i = d\-^ and 
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D 22 = d\ 2 represent intranode scattering at nodes 1 and 
2 respectively; D\ 2 = D 2 \ represent internode scattering. 

After projecting to the LLL, the disorder matrix re- 
duces to a 2 x 2 matrix and can be represented by the 
Pauli matrices rj. The most general disorder matrix pro- 
jected to LLL then takes the form 

^llL = (!0) 

3=0 

where we have denoted the / matrix by rjo ■ 



IV. AVERAGE DENSITY OF STATES 

Using a four component bosonic spinor <j) and a four 
component Grassmann spinor x the average retarded 
Green function for a noninteracting problem can be writ- 
ten as 

_ 3 

G R (E;r,r') = V[cj>*\D[<f\V[ X *]V[ X ] 

3=0 J 

V[V J ]P[V 3 W(r) ( j ) (r')e sR 1 (11) 
where P[Vj] is the probability distribution of Vj and 



= ijd 2 r 



b\E - H - H imp + iS)(j> + 



imp 



iS)x 



The average density of states is given by 



p(E) = --Im G R (E;r,r). 

7T 



(12) 



(13) 



After performing the disorder averages we can write 

G R (E;r,r') = -i j V[4>*]V[4>]V[x*)V[ X ] 

cj>*(r)cf>(r')e AR , (14) 

where the action A R involves interactions among the 
holds generated by the disorder averaging procedure. Af- 
ter projection to the LLL, the action A R can be expressed 
in terms of a two-component holomorphic bosonic spinor 



«'>-($))■ 



(15) 



and a two component holomorphic Grassmann spinor 

'wi{z) 



X{z) 



w 2 {z)l • 



(16) 



In terms of these holds the action is given by 



A 1 

A 
A 



A 



3=0 



f-I 



d ze 



d zhj 



-zz/2l\ 



-zz/11 



where e = E + iS and 



(17) 



(18) 



is the effective interaction of the helds generated by the 
averaging over the random variable Vj. For the Cauchy 
distribution, defined by 



9j 



1 



v 91 



we have 



h j( K ) = -9j\k\- 



(19) 



(20) 



If the disorder distribution is Gaussian white noise, de- 
fined by 



we get 



P[ Vj (f)}=AfeM-^- j J d 2 rV 3 2 (r}], 



(21) 



(22) 



The above action is invariant under the translation fol- 
lowed by a gauge transformation. Due to this invariance, 
the spatial dependence of the average retarded Green 
function is same as the spatial dependence of the pure 
system's Green function 



Gp U re( E > z li z 2) 



exp[ 



N 2 + N 2 



2z x z 2 ] 



2nl 2 B (E + iS) 



(23) 



The disorder averaged Green function can be written as 

G R (E, z u z 2 ) = C(E + iS, 9j ) exphdzi | 2 + \z 2 \ 2 - 2 Zl z 2 ], 

(24) 

where gj 's are coupling constants of various types of dis- 
order and C(E + iS,gj) is a gauge invariant proportion- 
ality constant which depends on the energy and disorder 
strengths. This gauge invariant proportionality constant 
is what we need to calculate to hnd the average density 
of states. 

For the calculation of the average Green function's de- 
pendence on the energy and disorder coupling constants 
we introduce two new Grassmann variables 9 and 9 and 
enlarge the Euclidean coordinate space into a superspace 
of coordinates (x,y,9,9). Integrals over the Grassmann 
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coordinates are normalized as ir J d9d999 = 1 . The norm 
of a coordinate vector is defined as x 2 +y 2 +99. This norm 
is invariant under the superspace rotations. In addition 
to the ordinary rotations in the Euclidean subspace and 
the symplectic transformations in the Grassmann sub- 
space, the superspace rotations involve transformations 
which mix (x, y) and (0, 9) in the following manner: 



f+ 2i x m + 2i 2 ne 

+ 4(f 2 -r)fi 
0-4(Z~i -f)n. 



(25) 



In the above set of transformations l\ 2 are two arbitrary 
Euclidean vectors and is a Grassmann number. We 
also define two holomorphic superfields and their conju- 
gates as 



$(z,9) 



(*) + 



${z,0) = <fi(z) + 



V2h 



V2h 



(26) 



In terms of these superfields the pure part of the action 
can be expressed as 



Af = 2iml 2 B / d 2 zd9d9e~^ +e ^/ 21 ^^, 



(27) 



which is manifestly invariant under superspace rotations. 
After the disorder contributions to the action are ex- 
pressed in terms of these new superfields, we have to 
demonstrate these to be invariant under superspace ro- 
tations. In order to be supersymmetric A^'s have to 

I 



be local in the supercoordinate space and this is only 
possible if they do not involve any quartic fermionic in- 
teractions. We note that 



e -zz/2? B ^ ri . (j) + x \ v . x ^ 



e zS/21b x^VjX 



e~ z ~ z/lB (x%x) 2 , (28) 



where hj and hj correspond to the first and second 
derivatives of hj with respect to its argument. The Tay- 
lor series truncates at the quadratic order as the higher 
powers of X VjX are identically zero according to the 
anticommutation rules. We also note that (x VjX) 2 — 
—2w\wiwlw2 for j = 1,2,3 and (x^Vox) 2 = 2w\w\w\w2- 
If hj does not vanish we get four-fermion interactions. 

If there were one bosonic and one Grassmann fields in- 
stead of spinors, as in the problem solved by of Brezin et 
al.r^ no four fermionic terms would be generated, and 
the action for an arbitrary disorder distribution would 
be local in the superspace coordinates. For the case un- 
der consideration, such a simplification is not possible in 
general. However, for Cauchy distribution the disorder 
averaged action is quadratic and can be made manifestly 
supersymmetric. Thus, the calculation of the DOS re- 
duces to a calculation of a zero dimensional field theory 
over two complex bosonic fields. 



A. Cauchy Distribution 



m — 



The action is given by 



A 1 



d 2 ze- zS ' 2l2 £ 



le 



+ x^x) -^2gj\^V3<t> + x ] mx\ 



(29) 



Using the superfields $ and $ the action can be written as, 



A* = 2nl 2 B I d 2 zd9d9e- (zg+eS V 21 * 



ie$$-J29j\®Vj®\ 

j=0 



(30) 



which is manifestly invariant under rotation and magnetic translation in superspace. Because of this symmetry, the 
DOS can be reduced to a simple expression involving integrals over two ordinary complex variables. Expressed in 
terms of two radial and two angular variables, it is 



~ p{E) = 2^K Im l ln 



roc p2tt p2tv 

d(r 2 /2) i d{rl/2) I dot! / da 2 exp 
it Jo Jo Jo 



i<r\+rl)- 9o\rl+r 2 2 \ 



-gzVi - rfl - 2gir 1 r 2 \ cos(«i - a 2 )| - 2g2r 1 r 2 \suv(a 1 - a 2 )\ 



(31) 
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For simplicity, consider the cases where we keep only one of the internode scattering, or the random mass term, 
along with the potential disorder. We get: 



(i) 9i =92 = 



(") 93= 92 = 0; 



p(E) = 



2n 2 l 2 B 



,9o 



,9o + .93 



.9o 2 + £ 2 



P(E) 



2n 2 l 



272 



(,9o+.93) 2 + £ 2 
.90 + .9i 



go 

g 2 + E 2 + (g Q+gi f + E 2 



(32) 



(33) 



The answer for the case (Hi) g$ = g\ = is identical to the case (ii). The DOS obtained for these three cases are 
identical, as the Hamiltonian involves only one Pauli matrix at a time, and these matrices are related by unitary 
transformations. 

Consider now g\ — g 2 = giN and 53 = 0. We obtain, defining by / the expression within the curly parenthesis in 
Eq.EB 



I =- 



(« 2 - Hn) 



9in 



9 2 in 



: tan 



9IN 



V^ 2 



9 2 in 



2\/27r— + 2tt 



9in 



9 2 in 



(34) 



where a = g — ie. The details of the evaluation of the multiple integrals are provided in the Appendix [X] The 
expression for the DOS obtained from this expression is lengthy and not very illuminating, but it is important to note 

that because of the presence of the term tan -1 (c/in l^/ot 2 — djpfj > we obtain a ln_E divergence at the band center 
when go = 0- Based on symmetry, similar behavior will be obtained when a combination of two Pauli matrices are 
considered. This should be contrasted with the (InE) 2 divergence obtained by Hikami et 



■ 55 



2. m# 



When the fermion is massive, we will ignore the mass disorder part. The density of states is given by 



P(E) 



1 



Im 



_8_ _d_ 

dei de 2 



In 



ci(r 2 /2) / d(r|/2) / da x / da 2 exp 



l r i + r l\~ 2gir-ir 2 | cos(ai - a 2 )\ - 2g 2 r 1 r 2 \ sin(ai - a 2 ) 



(35) 



where ei 2 = E ± m + iS. Again if we take only one of the internode scattering terms (g 2 = 0) for simplicity, the 
expression within the curly parenthesis in Eq. 1351 /, becomes 



The density of states is then given by 



An 2 1 



2/2 ^2 



.90 



ab + gwfab 

1 g (Rcosf3 + giVficosf ) + E(Rsinf3 + c/iV^Rsinf ) 



R 2 + g 2 R + 2 Sl Ri cos f 



(36) 



(37) 



where R = y/ (g^ + m 2 — E 2 ) 2 + Ag 2 E 2 and tan/3 = (2g^E/ (g 2 , + m 2 — E 2 )). The above expression takes particularly 
simple form when go = 0. It becomes 



p(E 2 > to 2 ) 



1 



47^1 L 



5(E + m) + 5(E - to) + 



2E 9l 



(E 2 -m 2 + g^VE 2 



p(E 2 < to 2 ) = 0. 

If both internode scatterings are present and g\ — g 2 — giN, the integral is given by 



(ab - 2gj N ) 



n - 4 



9in 



tan 



1 I 9in 



\/ab-g 2 IN \\/ab-g 2 IN 



2V2n^ + 2n- 9lN 



ab 



y/ab-g 2 IN 



(38) 



(39) 



The expression for the DOS is tedious. However for g — 0,the feature that the DOS is zero for E 2 < to 2 is still valid. 
In this case for energies close to ±to, p(E) ~ In \E — m\/ \J\E — m\. 
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V. HALL PLATEAU IN THE LOWEST LANDAU LEVEL 



Similar to the method described in Ref. l58l . we generate the matrix elements of the Dirac Hamiltonian after 
projecting to the lowest Landau level. In our problem, the element (k\H\k') itself is a 2 x 2 matrix: 



\H\k') 



dxdyipl (x, y) [H 1 ^ (x,y)+ mi] 3 ]ip k/ (x, y) = mr) 3 5 ktk/ + V(k, k') 



(40) 



r 



where ip k (x,y) is the lowest Landau level wave function elements of the 2x2 matrix V(k, k') can be computed 

in the Landau gauge. We choose all the Vj's to follow explicitly — for example, 
independent Gaussian white noise distributions such that 
V 3 {x,y)V r {x',y') = g^5 j .,&{x - x')5{y - y'). Then the 



V{k,k') lx = 



^ \^ B (fc ) 



J d£ e ^{gouoiht 



^l%,k> -k) + g 3 u 3 {l B i + ^-^1%, k' k)], (41) 



where Uj(x, k) is a complex random variable defined to It is straightforward to compute the statistical proper- 
be the Fourier transforms of Vj (x, y) along the y direction ties of the matrix elements. The averages are: 
normalized by the width gj, namely: 



Uj (x, k) 



9jV L 



dyV 3 {x,y)e lk y. (42) 



Because each of the disorder fields has zero correlation 
length, and there are no correlations between them, 



HA'. A'':,.. : i,j = 1,2 



(44) 



Ui(x, k)uj(x' , k') = SijS(x — x')5(k + k') (43) As to correlations, the only non-vanishing pairs are: 
I 



V{k 1 ,k 2 )uV{k 3 ,k i ) 11 


= V(ki,k 2 ) 22 V[k 3 , fc 4 ) 22 = 


V{kl,k 2 ) n V{k 3 ,k 4 ) 22 


= V r (fcl,fa)22V r (& 3 ,fc4)n = 


V{k l ,k 2 )i2V{k 3 ,k i )xi 


= V(k 1 ,k 2 ) 21 V{k 3 ,k i ) 21 = 


V{k 1 ,k 2 )l2V{k 3 ,k i ) 2 l 


= V(k 1 ,k 2 ) 21 V(k 3 ,k i ) 12 = 



9o + .93 



2ttL, 

1 1 

% - .93 



exp | 



2 2 

91 - 9 2 



2nLy 
9 2 i + 9l 



2nL,, 



exp | 



exp | 



I 2 


-k 2 y~ 


- (ki 


-ki) 2 )]S kl 


I 2 


-k 2 ?~ 


- (fc 4 


-ki) 2 )]S kl 


l 2 


-k 2 f- 


- (fc 4 


-ki) 2 )]5 kl 


l 2 


-k 2 f- 


- (fc 4 


-ki) 2 )]6 kl 



(45) 



For numerical implementation, we discretize and use <5-correlated as in (|43[) . Finally, we approximate the in- 
the integer I to label the x coordinate. We then gener- tegrals by sums. Explicitly, the matrix elements are 
ate a set of complex random variables Uj(l,k), that are 
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V(k,k + k') n 
V(k,k + k') 22 
V(k,k + k') 12 
V{k,k + k') 21 



2 , 12 

-a k 



2 i '2 

-ark 



2 , 12 

~—a k 



,— a fe 



VMA 



J2 [9oUo(2k + j, k') + g 3 u 3 (2k + j, k')} e- a2 > 2 

3 

[9ou (2k + j, k') - g 3 u 3 (2k + j, k')] e~ a ^ 2 

3 

^2 [ 5 iMi(2fc + j, k') - ig 2 u 2 (2k + j, k')] e- a ^ 2 

3 

[ffiui(2* + j, k') + ig 2 u 2 (2k + j, k')] 



(46) 



with A = V . e~ 2a 3 and a 2 = n/2M 2 . Here M is the duce fluctuations. Energies {Ei} were chosen close to the 

critical energy and measured in units of 2(^ J ._ p?) 1 / 2 



length of the system in the y direction, the unit being 
y/2nlB, that is, M = L y / ^/2ttIb, chosen to be an inte- 
ger. The integers k and k' label the wave vectors. Since 
the matrix elements decay exponentially, we can neglect 
them for k 1 > 2M . A cutoff is also necessary for the 
recursive Green's function technique— that we use. 

We compute the density of states p(E) by directly diag- 
onalizing the Hamiltonian. We have checked that p(E) 
is independent of M, for sufficiently large M; M = 32 
seems to be sufficient; the total number of momentum 
states Nk is chosen to be 1000, which is half the dimen- 
sion of the Hamiltonian matrix to be diagonalized, as 
there are two fermions for each k. Typically, an average 
over 100 disorder realizations is used. 

The recursive Green's function technique, similar to 
that in Ref. l58l . is used to explore the localization prop- 
erties. The details are described in Appendix [B] We 
first compute the localization lengths for a finite sys- 
tem, XM(Ei), at a set of energies, in systems 



with transverse dimensions {My} Since there are two 
types of fermions, in general there can be two distinct lo- 
calization lengths; however, in most cases discussed be- 
low, they are identical within our numerical accuracy, 
and we will not generally distinguish them. Assuming 
finite-size scaling, \ M (E)/M = f(M 1 /' y (E-E c )), where 
fix) is a universal function, the data is collapsed to ob- 
tain the localization length exponent is, and the critical 
energy E c . Strictly, scaling holds only for large enough 
systems in the vicinity of critical energy. Here the ener- 
gies {Ei} are chosen close to the critical energy, E c , and 
the validity of the scaling law is verified by the success of 
data collapse. For the details of the procedure involving 
data collapse, see Appendix ICl 

The numerical calculations about the localization 
properties were mostly performed for a quasi-one di- 
mensional system with the transverse dimensions M — 
8, 16, 32, 64. The total number of momentum states is 
Nk = 5 x 10 4 . Because of the 2x2 character of the 
Hamiltonian matrix elements, the numerical calculations 
are more demanding than those in Ref. HH. The data are 
typically averaged over 100 disorder configurations to re- 



similar to Ref. |5£ 

Our program is also validated by the case go — 0.5, 
9i = g 2 = g 3 = 0, and m = 0. In this case, the two types 
of fermions are independent. Because the LLL wave func- 
tion is identical to the non-relativistic one, the properties 
should be the same as in Ref. l58l . Numerical computa- 
tions show a single peak in the density of states and a 
localization length exponent of v = 2.41 ± 0.08; both 
agree well with the previous results. 



VI. LD TRANSITION FOR THE MASSLESS 
CASE 

A. One disorder field 

Consider first the cases where only one type of disorder 
has nonzero strength. Numerically, we considered (1) 
.93 = 0.5, ,g = gi = g 2 = 0, (2) g 1 = 0.5, g = g 2 = g 3 = 
0, and (3)^2 = 0.5, go = g\ = g 3 = 0. In all of these 
cases, the delta function density of states in the pure 
system are broadened into a simple bell shape function 
due to disorder. The results of successful data collapse, 
not shown here, yield critical exponents v = 2.46 ± 0.09, 
v = 2.48 ± 0.11 and v = 2.45 ± 0.08 for cases (1), (2), 
and (3), respectively, that is, they are the same within 
the error bars. 

The critical exponents are all equal to that of sin- 
gle type of fermions subject to potential disorder. This 
can be understood as follows. The original Hamiltonian 
matrix is in the basis {|/ci, 1), |fci, 2), \k 2 , 1), \k 2 , 2), . . .}, 
where 1 and 2 label the type of fermions. If we re- 
order the basis as {|&i, 1), \k 2 , 1), . . . , \ki, 2), \k 2 , 2), . . .}, 
the Hamiltonian becomes a 2 x 2 block matrix with di- 
agonal blocks representing intranode parts, and the off- 
diagonal blocks representing internode scatterings. Ex- 
plicitly, it is in the form: 



H = 



9qUq + g s U 3 + ml g\U\-ig 2 U 2 
giU!+ig 2 U 2 g U - g 3 U 3 - ml 



(47) 
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where [/,-, i — 0, ... ,3 are statistically independent ran- 
dom Landau matrices, and I is the identity matrix. The 
case with only g% nonzero, has the same structure, and 
hence same statistical properties, as the case when only 
go is nonzero. When only g\ is nonzero, we can, by a 
unitary transformation given by, 



T = B 



I I 

I -I 



(48) 



where B is a normalization constant, bring the Hamilto- 
nian back to the block diagonal form, resulting in a struc- 
ture corresponding to two types of independent fermions 
in the presence of mass disorder. Thus, the critical ex- 
ponent is the same as the case when only g is non-zero. 
The same argument also applies to the case when only 
(72 is non-zero. 



B. Two disorder fields 

If one of the two disorder fields is Vo(x, y), and another 
is V±(x,y), or V2{x,y), or V^(x 1 y), an appropriate uni- 
tary transformation about an axis by tt/2 will map one 
possible case to another. For instance, the transforma- 
tion in Eq. (|48p will transform the case with Vo(x, y) and 
V2(x,y) to Vo(x,y) and Vs(x,y). It is therefore sufficient 
to consider only the case with just Vo(x, y) and V^{x,y). 
However, from Eq. (|47|) . the Hamiltonian is block diag- 
onal, and the blocks goUo + 53^3 and goUo — 93U3 are 
statistically equivalent to a new block -^/ffo + diU, with 
U a new random matrix satisfying the same statistical 
properties as U^s; see Eq. (|45|) . That is, the Hamiltonian 
for 50 7^ and g% 7^ is statistically the s ame as t hat 
corresponding to a potential disorder g = \/ g$ + g\ ■ 

Our numerical computations confirm this argument. 
The data collapse was found to be successful, assuming 
E c = 0, and the critical exponents are v = 2.45 ± 0.06 
for g = 53 = 0.5, gi = g2 = 0, and v = 2.45 ±0.11 for 
9o = 9i = 0.5, g 2 = 33 = 0. 

Next, we choose two disorder fields from V\{x,y) 1 
V2(x,y), and V^{x,y). There are three possible combina- 
tions. In fact, these three cases are not independent; we 
can map one case to another by an appropriate unitary 
transformation corresponding to a rotation by tt/2 about 
a certain axis. Therefore, it is sufficient to consider only 
one of the three cases; for example, let us choose Vs(x, y) 
(mass disorder), and V\{x,y) (internode coupling). 

We set gi = gz = 0.5, go = 92 = 0. The density of 
states and the localization lengths are plotted in Fig. [1] 
for E > 0; there is symmetry under E — > — E. The 
extended states are no longer at E — but shifted to 
E = E c ~ ±0.42. At E c ~ ±0.42, we study the localiza- 
tion properties using the data in the range of \E\ > E Cl 
since data in the range \E\ < E c are close to both critical 
points and are likely to result in inaccurate results. The 
maximum system size used is M = 64. Because we do 
not have a priori knowledge of E c , the statistical proce- 
dure discussed in Appendix [C] is employed to determine 




FIG. 1: (Color online) The DOS (top) and the localization 
length in finite systems (bottom) for gi = gs = 0.5, go ~ g2 = 
0. The dashed line shows the LD transition. 



1 ■ «iiT 

* 7 H 



■ M=8 

• M=16 

* M=32 
T M=64 




-M=8 
-M=16 
-M=32 
-M=64 



M (E-E c ) 

FIG. 2: (Color online) Scaling curve for the case gi = g-j, — 
0.5, and go = g2 ~ 0. Insert: the dependence of Xm/M on 
the energy E for different system sizes M. The shaded area 
is used for scaling and the critical exponent is found to be 
v = 3.23 ±0.26 



E c , hence the critical exponent v. The data collapse is 
shown in Fig. O The critical exponent for this parame- 
ter set is found to be v — 3.23 ± 0.26, distinct from the 
nonrelativistic case of v ~ 7/3. 

The present problem can be exactly mapped onto the 
the spin-orbit scattering involving the two-state Landau 
level problem discussed in Ref. [6l|; our results are in full 
agreement. From Fig. [TJ there appears to be a diver- 
gence in the DOS at the band center, corresponding to 
a possible LD transition at E — 0. As shown above, the 
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0.8- 




-1.5 -1.0 -0.5 0.0 0.5 1.0 1.5 

E 



FIG. 3: The DOS when the mass disorder gz and the in- 
ternode coupling gi are both present. Parameters are: top: 
gz — 10gi = 0.15; middle: 33 = gi = 0.15; bottom: 
53 = O.lffi = 0.15. 



Cauchy distribution does lead to a In E divergence, but 
such a weak divergence is difficult to detect numerically; 
note, however, that the numerical calculation involves 
Gaussian disorder. A semiclassical explanation^ of the 
existence of an extended state at the band center for the 
two-state Landau level problem is known. However, such 
an argument is delicate and fails if a the third kind of 
disorder is present, which is likely in graphene, where 
potential disorder can not be avoided. Thus, we shall 
not consider further the possible extended state at the 
band center. 

It is interesting to study the behavior as the ratio g\ / g% 
is varied. The result for the DOS is show n in Fig . [3] 
Note that the energy E is in the unit of 2 \J g\ + g\ ■ So 
the extent of the band increases, as g\ increases. The 
divergence of the DOS at the band center is a unique 
feature when 171 and (73 are equal, while in the extreme 
limits there may be a slight dip at E = 0. 

As to v, a continuously exponent is suggested in Fig.|U 
In the limits 173 g\ or 53 -C g\ , only one type of disor- 
der dominates, hence the value v ~ 7/3 is plausible. The 
deviation from this value is the largest when gi ~ 173, 
although the data collapse becomes insensitive to the 
value of v in the same regime, resulting in larger error. 
Nonetheless, the results are suggestive of a continuously 
varying critical exponent. 

In the Hamiltonian, the mass disorder V^{x, y) and the 
internode scattering disorder Vi(x,y) are accompanied 




FIG. 4: The dependence of the exponent v on gi/gz- The 
parameters are go — g^ = 0, gz = 0.15, and m = 0. The 
dashed line corresponds to v — 7/3. 



by the Pauli matrices 773 and 771. If we apply a unitary 
transformation 

T = c(j- J ) (49) 

corresponding to a rotation of tt/2 about y axis, where C 
is a normalization factor, the disorder Hamiltonian (147)) 
will be transformed such that (73 — > gi, g\ — > —(73. Be- 
cause we are studying statistical properties of the system, 
and all distribution functions are symmetric about zero, 
the negative sign in front of the 33 is of no importance. 
This means that this unitary transformation effectively 
interchanges g\ and 53 , hence map the regime 173 > g\ to 
the regime 173 < g\. Note the symmetry between the two 
regimes in Fig. [3] and Fig. [4j 



C. Three disorder fields 

The important case in this category is when g\ ~ g 2 ~ 
173; other cases can be roughly understood in terms of 
the cases discussed above. For numerical computation, 
we take g\ — g 2 = gz = 0.5. The DOS is shown in Fig.O 
Compared to the case when only g\ = gz — 0.5, discussed 
above, the divergence of the DOS at E = is missing, but 
the two peaks at E = ±0.46 survive. This is suggestive 
of nonexistence of extended states at the band center, 
but a LD transition at E ~ ±0.46, which is confirmed by 
the scaling curve shown in Fig. [HI and a critical exponent 
of v = 3.6 ± 0.3 is obtained. The error bar is large due 
to a substantial degree of disorder, but the exponent is 
distinctly different from the value v = 7/3, indicating a 
new universality class. 
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FIG. 5: The DOS as a function of energy E in the case g\ = 
92= g-i = 0.5, go = 0. 





0.5 1.0 



0.1 1 

M""(E-E c ) 

FIG. 6: (Color online) The scaling curve for the case go = 0, 
and gi = g2 = ga = 0.5. Insert: the dependence of Xm/M 
on the energy E for different system sizes M. The shaded 
area is used for scaling. The critical exponent v is found to 
be 3.6 ±0.3 



D. Four disorder fields 

Finally, we have also examined the case go ^ 0, g\ ^ 0. 
32 7^ 0, <?3 ^ 0. We have found that when g is small, the 
potential disorder simply broadens the density of states, 
and results in a value of v as though it did not exist; on 
the other hand a large value of go drives v to a value close 
to 7/3. 



VII. LD TRANSITION FOR THE MASSIVE 
CASE 

The constant mass term results in new physics when 
combined with the internode coupling, and the resulting 
phenomena are different from the case when the mass 
disorder and the internode coupling are combined, as in 
the previous section. 




0.4 0.6 0.8 1.0 1.2 1.4 1.6 1.8 2.0 2.2 



FIG. 7: The DOS as a function of energy E. In the case go = 
<?2 = 33 = 0, g\ = 0.5, and m = 0.5, the DOS vanishes when 
\E\ < m. Insert: the logarithmic plot of the DOS around E = 
m, and m = 0.5. The best fit gives a slope of —0.523 ± .003. 



A. One disorder field 

Consider first the case m ^ 0, and either go or 173 is 
nonzero. From Eq. (]4T[) . the random Hamiltonian matrix 
are block diagonal, hence the two types of nodal fermions 
are uncoupled. The energy of the two fermions are shifted 
by the amount of ±m, and the DOS is simply a super- 
position of two bell shaped functions centered at ±m. 
As to localization properties, because the two fermions 
will have different critical energy, namely E c = ±m, at a 
particular energy E these localization lengths will be dif- 
ferent from each other. The data collapse for the fermion 
with E c — m once again gives a critical exponent of 
2/~7/3. 

Consider now finite mass m ^ and only one intern- 
ode coupling, for example g\ 7^ 0. The calculated DOS 
with E > region is shown in Fig. [7] for the parameter 
set go = g-2 = .93 = 0, gi = 0.5, and m — 0.5. Qual- 
itatively the results are similar to the analytical calcu- 
lations involving the Cauchy distribution, even though 
the numerical computation is for the Gaussian distri- 
bution of disorder. First, the DOS vanishes in the re- 
gion \E\ < m. Second, the analytical calculation shows 
that p(E) ~ (E — m) -1 / 2 , E — > m + , independent of the 
value of g\ . The insert of Fig. [7] yields an exponent of 
— 0.523 ± .003. We have checked that this value does not 
vary with gi, within our numerical accuracy. 

As to LD transition, we perform data collapse with 
E c = m. The critical exponent turns out to be v = 
1.82±0.06 for the parameters go = 52 = .93 = 0, .91 = 0.4, 
and m = 0.15, see Fig. [SI which is significantly differ- 
ent from the usual case corresponding to v ps 7/3. This 
striking result implies that the system belongs to a new 
universality class. We now vary g\, keeping m fixed to 
0.15, and the result for the exponent v is shown in Fig. [SI 
It appears that the exponents continuously vary with the 
ratio gi/m. 
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FIG. 8: (Color online) The scaling curve for the case go = 
52 = 33 = 0, gi = 0.4, and m = 0.15. Insert: dependence 
of Xm/M on energy E for different system sizes M. The 
shaded area is used for scaling. Critical exponent is found to 
be v = 1.82 ± 0.06, with the choice of E c — to. 




FIG. 9: The dependence of the critical exponent v on gi 
(normalized by the mass to). The parameters are go = <?2 = 
gz = 0, and to = 0.15. The dashed line indicates the level of 
v = 2.33. 



B. Two disorder fields 

The most relevant case corresponding to graphene is 
the one with two types of internode scattering of com- 
parable magnitude. Therefore, we choose <?i = <?2 = 
to = 0.5. The DOS is shown in Fig. [TOl Note t hat the 
energy is now measured in units of 2 \J g\ + g\ = V2, 
so that the divergence is located at E = m, which is 
m = 0.5/V2 = 0.354. From the insert in Fig. [TU1 we 
find a slope of —0.47 ± .01, which is consistent with 
the analytical result for the Cauchy distribution, namely, 
p(E) ~ In \E — m\j ^J\E — m\. Also note the gap in the 
DOS. 

The crossing point in Fig. [10] indicates E c ~ 0.55 in- 
stead of E c = to. The data collapse is shown in Fig. [11] 
with a critical exponent of v = 3.8 ± 0.2. This crit- 
ical exponent is not close to any of the values found 




FIG. 10: (Color online) The DOS (top) and the localization 
length in finite systems (bottom) for g\ — g^ = m = 0.5, 
9o — <?3 = 0. Insert in the top: the logarithmic plot of the 
DOS around E = to. The best fit gives a slope of — 0.47± .01. 



for a finite mass with a single internode coupling, as in 
Fig. [HI However, it is reasonably close to the exponent 
for to = go = gi = 0, and g\ ~ 33 (see Fig. 2]), which is 
also equivalent to the case m = go = gs = 0, and g\ ~ g?, 
as discussed above. No te that to = 0.5 is much smaller 
than the bandwidth 2y g\ + g\ = 1.414, and this critical 
exponent indicates that the presence of small finite mass 
will have little effect on the critical exponent as long as 
two internode couplings are finite. 

At the other limit, when m is large enough compared 
to the band width, the critical exponent v — * 1. Thus, 
it is reasonable to believe that the exponent also varies 
continuously, as a function of <?i/to, provided that g\ — 
</2, and the behavior is similar to in Fig. [8] except that 
v — > 3.8 in the limit to — > 0. 



C. Three disorder fields 

Because potential disorder is always present in exper- 
iments on graphene, we would like to discuss the case 
when go 7^ 0, to ^ 0, and g\ = gi 7^ 0. 

First, when to is the smallest parameter, it can be 
neglected, and hence the massless case discussed above 
is recovered. Our numerical computations gave a critical 
exponent of v ~ 3.8 when to <§C go <§C <?i = <?2, and 
v ~ 2.3 in the limit to <C <?i = 32 *C go, because now the 
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FIG. 11: (Color online) The scaling curve for the case gi — 
g2 — m — 0.5, go = gs = 0. Insert: the dependence of Xm/M 
on the energy E for different system sizes M. The shaded 
area is used for scaling. 

potential disorder is more important than the rest. 

When go is the smallest, it will have little effect. There- 
fore, as discussed in the previous subsection, there will 
be a continuously varying exponent from v ~ 1.0 for 
go < g\ = gi < m to v ~ 3.8 for g < m -C gi = 32- 

Finally, if 171 = 32 are smaller than the rest, the in- 
ternode scattering is no longer important, and hence the 
two nodal fermions will be decoupled. The exponent will 
therefore be always v ~ 2.3 regardless of the relationship 
between m and go. 



ternode scattering strength is larger than the intranode 
scattering strength, we expect that in the lowest Landau 
level — > ±1, 1 — > 2 and —1^—2 plateau transitions 
can have different universality classes, in contrast to the 
conventional IQH effect. In the opposite limit, when the 
intranode scattering is considerably stronger, the plateau 
transitions will fall into the conventional IQH universal- 
ity class with v ~ 2.3. 

For the spin and nodally degenerate plateaus, the po- 
tential scattering is strong compared to the Zeeman en- 
ergy and the mass gap. Theoretically, from our analysis 
of the massless cases we can infer the plateau transitions 
to be of the conventional IQH type. However, in exper- 
iments, if scaling with respect to the band center is in- 
voked, an effective exponent for these plateau transitions 
will be observed. 

Plateaus at Vf = ±4, ±6, .. involve higher Landau lev- 
els. In the higher Landau levels both the random poten- 
tial at the lattice scale and the random hopping will have 
nonzero intranode as well as internode scattering contri- 
butions. This will complicate the analysis of LD transi- 
tions in these levels. Because inclusion of a finite mass 
does not lift the nodal degeneracies of higher Landau lev- 
els, the effect of internode scattering can be strong. In 
view of the degeneracy factor, there is the possibility of 
observing an effective exponent for these higher plateau 
transitions. 



VIII. CONCLUSIONS 



We have analyzed the effects of disorder on the LD 
transition in the LLL of graphene. Because both types 
of internode scattering, present in the LLL, arise from 
the random hopping, they will have roughly the same 
strength. Because the sources of the mass disorder and 
the internode scattering are different, their strengths will 
be generically different. In some special cases of disorder 
combinations we have found new universality classes of 
LD transition in contrast to the conventional IQH. 

Our results for the LD transitions in the LLL have di- 
rect experimental relevance for the plateau transitions in 
graphene. Consider first the cases where both the spin 
and the nodal degeneracies are completely removed. A 
number of authors have shown that the inclusion of a fi- 
nite mass and Zeeman energy can explain the appearance 
of plateaus at Uf = 0, ±1, ±2q. 

Because experiments resolve the spin and nodal split- 
ting, intranode scattering which always broadens the 
Landau levels is weak compared to E z and m. If the in- 
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APPENDIX A: INTEGRAL 

For the Cauchy distribution, in the absence of mass 
disorder, the calculation of the DOS in Eq. (|35|) involves 
the integral 
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I dx 
'o JO 



dy 



2tt p2ir 

da\ I da.2 exp 
it Jo 



by - 2g 2 ^/xy\ sin(ai - a 2 )| - 2g 1 ^/xy\ cos(«i - a 2 )| 



(Al) 



r 



where a = g$ — ie\ and b = go — ie 2 . In the massless 
case a — b — go — *£• After performing the integrals over 
one of the angles,the double integral over the angles is 
reduced to 



L ang 



/■t/2 

8vr / daiCxp[-2A7iycos(ai - /?)], (A2) 
Jo 



where A = y^f + gf and tan/3 = 52/31 ■ Now, expand- 
ing the exponential in a power series, integrals over x and 
y can be easily performed. For the angular integral we 
use the relation 



tt/2 



da 1 cos'(ai — (3) 



1+1 



i ;+1 M(^,i^,sin 2 /3) 



l+l a v f l + 1 1 3 + ^ 2 m 
cos' +1 ^2-^i(— 1 2' ""a - ' cos $ 



(A3) 



where 2-F1 is Gauss's hypergeometric function and obtain 
Aa6 



J = 



2;r ^^ r 2 (J/2) + 

ra + 2) 

j=l ;=o v ; 



1) ."2<?w 



_ J + l 1 3 + 1 2 2 

x 2 i J, i(—, 5,— ,#/A 2 ). 



(A4) 



The summation over Z can be performed by the follow- 
ing trick: (1) use the integral representation of 2-F1, (2) 
perform a power series summation which is simple in 
these cases, and (3) complete the integration over the 
auxiliary variable introduced for the integral representa- 
tion. For simplicity we will specialize to the cases: (i) 
.9i 7^ 0, 32 = and (ii) gi = g 2 = giN- 
(i) <7i 7^ 0, g 2 = 0: In this case we have 



r 2 ((Z/2) + 1) -2ffi 



;=o 



r(i 



(-f=) 



j + i 1 3 + ; 

x 2 ^(— ,1) 

We now use the following two relations 
1 

2 '2' 2 



J + l 1 3 + Z v^r((3 + z)/2) 

2-^ 1 V o ' o> ~~ ' i ' _ 



z=o 



r((z/2) 

< r((f/2) + i)r((3 + Q/2) 
r(z + 2) 



2 + x 



(A5) 

(A6) 
(A7) 



to obtain Eq. 

(ii) 31 = 32 = giN- In this case we have 



i=0 



r 2 ((z/2) + i) ,-2 5/JV ,, 



r(i- 



-)' 



y 2 '2 2 77 



(A8) 



Using the integral representation 



_ J+l 1 3 + / 1 

2^l(— JJ— , ~, — , 



£((3 + 0/2) 
F(l/2)r((//2) + l) 7 

x (^_.)V2 (1 _i)-l/2 

V l-</2 y V 2^ 



(A9) 



and Eq. (|A7|) we get 



J = -irVzVab [ 
Jo 



dt 



faby/t(l - (t/2)) + ,9/at Vt(l-i) 

(A10) 

After performing the integral over t we get Eq. (1551) . Af- 
ter setting m = one recovers Eq. 



APPENDIX B: RECURSIVE GREEN'S 
FUNCTION 



Because there are two types of fermions in our problem, 
the recursive Green's function technique is a little more 
complicated than that introduced by Huckenstein^. All 
the matrix elements, such as those in Eqn. (6.9) in 
Ref. [H, become 2x2 matrices; hence, all the operations, 
such as the multiplication and inversion are matrix oper- 
ations. 

Denote the 2x2 matrix (i\G\j) simply by G(i,j). Sup- 
pose G( K )(i,j), i,j = 1,...,K, the Green's function 
containing K momentum states is available, then as we 
add another momentum state, the recursion relations for 
G (K+1) (i, 3) are: 
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G (K+1) {K + 1,K+ 1) = \E-V(K + l,K+l)-^2V(i,K +\)^G (K \i 1 i)V{j,K + 1) 



h3 



G< K+1 \i,K+l) 



Y,G( K \i,j)V(j,K + l) 



G {K+1) {K + 1,K+1) : i<K 



(Bl) 



G^ +1 ) (i, i) = G<*> + G^ +1 ' (i, A" + 1)G( K+1 ) (K + 1, K + l)- 1 ^ 1 ) (A + 1, j) : i,j < K 

These matrix inversions here can be accurately computed because the sizes are small. Corresponding to the two 
types of fermions, we are interested in G"(l, K) nn , n = 1,2, because the localization length of a system of width 
M for the n th fermion, Am n , is related to this quantity by, 



\-i _ _ 

A M,n — 



M 



M 



ln|GW(l,A)„, n | 



A" 



KV2ttI b 



fc=i 



(B2) 



where g4 = G^ K \l, K) n>n /G^ K 1 ^(1, A — l) n ,n- Furthermore, if we define a set of 2 x 2 matrices g^ K \j), j — 
1,..., A, such that their elements g^ K '(j) m , n = G^ K \\,j) m ^ n /G^ K >{\,K) m ^ m , from the recursion relations (|B1[) . we 
obtain: 



-,(K+1) 



9 {K+1) ii)r 



^9 {K Hj)V(j,K + l) G^ K+1 \K + 1,K+1) 



1 



;/ K 'V>„,„ [^^(jM^ + i)g (k+1) (^ + i, 



i = K + l 



i < A (B3) 
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APPENDIX C: DATA COLLAPSE 

In this appendix, we describe how we can extract the 
the exponent v, and the critical energy E Cl if necessary, 
based on the computed localization lengths in finite sys- 
tems Am (E) assuming a single parameter scaling assump- 
tion. 

Suppose we have obtained {Xm(E)/M} in systems 
with {Mi}^J{ for {Ej}^^, each with a standard devi- 
ation {cr^E.}. Our goal is to find out the proper values 
of v and E c such that all the Ne X Nm data points col- 
lapse on to a single curve: 



Xm(E) 
M 



f(MV v [E-E c 



(CI) 



Since f(x) is unknown, it is difficult to characterize the 
quality of the data collapse. To overcome this difficulty, 
we proceed as follows: suppose that we are given a pair 



of values (y,E c ), we can attempt to represent the un- 
known function f{x; v 1 E c ) by a polynomial of degree N 
by simply performing a general fit to Eq. (|C2[) given be- 
low, based on a total of Ne x Nm data points {(x, y,<r y ) : 

Qog[Mj /v {Ei - A c )],log[A M (^0/^iW^)} : 



log 



Am (E) 
M 



N 

]Ta fc [logiM^iE-Ec)) 



(C2) 



where {a,i}fL are the coefficients to be fitted. In the 
computer implementation, the order of polynomials was 
chosen to be N = 5, since no significant changes were 
noted by increasing N to 9. The quality of this fhp^ is 
represented by the variable S defined as: 



N E x N M 

S(u,E e )= £ 



V% - f(xi,v,E c ) 



1 2 



(C3) 



If the preset values (v, E c ) are not correct, the data points 
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will be scattered, resulting in a large value of S, which 
in turn indicates a poor data collapse. However, when 
(is, E c ) attain the correct localization length exponent 
and the correct critical energy, respectively, S will be 
minimized. Following this procedure, by minimizing S 
with the standard gradient descent method, we are able 
to determine correctly both the critical energy E c and the 
localization length exponent v. Because the scaling law is 
only valid in the close vicinity of the critical energy, once 
E c is obtained from above procedure, we have to check, 
for the purpose of self consistency, that all energies used 
in the data collapse are indeed close to E c . 

As for the statistical error of u, the usual procedure is 
to assume that the minimized S m in follows a x 2 distribu- 
tion, and hence the error bar can be drawn corresponding 
to a certain confidence probability. However, this is not 
the case in this problem, since S m in does not follow the 
X 2 distribution due to the nonlinear form of the estimated 
parameters v and E c in (|C2|) .^ To draw an error bar for 
v statistically correctly, recall that we have the original 



data {(xi,yi,a Vi )}. We generate a large number of data 



„(*0 „(*) Jk) 



a£')}, fork = 1,2, 



sets synthetically {(x\ ' , y- 

such that x$ — x% , a^y) = a Vi , and y\ k ' a variable ran- 
domly distributed in the Gaussian form with a mean of 
yi and a standard deviation of <j yi . Next, we perform 
exactly the same procedure to get for each synthetic 

data set {(x[ k \ yf\ as was performed in actual 
data set {(xi,yi,a Vi )} for estimating the v and E c . Fi- 
nally, the error bar for v is drawn as the estimated stan- 
dard deviation: 



N, 



1 Ns 



k=l 



1/2 



(C4) 



where N s is the number of synthetic data sets, and v^i 
is the average of v^ k \ N s = 10 4 in computer implemen- 
tation. 
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